{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "\n",
    "## <center>Открытый курс по машинному обучению. Сессия № 2</center> ##\n",
    "# <center>Индивидуальный проект по анализу данных </center>#\n",
    "### <center>Автор материала: Кирилл Власов (@vlasoff)</center> ###\n",
    "\n",
    "\n",
    "\n",
    "##  Часть 0. Импорт библиотек"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as st\n",
    "import sklearn\n",
    "from sklearn import linear_model,metrics,model_selection, preprocessing\n",
    "from tqdm import tqdm_notebook \n",
    "from matplotlib import pyplot as plt\n",
    "import geopy as geo\n",
    "import geopy.distance as dist\n",
    "import arcgis\n",
    "from arcgis.gis import GIS\n",
    "import folium\n",
    "from folium import plugins\n",
    "import plotly.plotly as py\n",
    "import plotly.graph_objs as go\n",
    "import seaborn as sns\n",
    "\n",
    "#Отключим предупредения\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 1. Описание набора данных и признаков"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Цели и задачи ###\n",
    "Спрогнозировать кол-во посещений в музеях РФ в 2016 году по историческим данным. Сделать выводы о влиянии различных факторов на посещаемость. \n",
    "> Прим.: Известны данные с 2008 по 2016 год, поэтому будем прогнозировать именно 2016, однако, \"отрежем\" данные за этот год, позже - после некоторых преобразований\n",
    "\n",
    "### Источник данных ###\n",
    "\n",
    "Для исследования использованны статистические данные по музеям с официального сайта <a href =\"http://opendata.mkrf.ru/opendata/7705851331-stat_museum/\">Министерства Культуры России</a>\n",
    "Данные за последние 9 лет можно скачать в удобном для работы csv формате, однако сайт не позволяет скачивать данные сразу за несколько лет, поэтому на входе получены 9 файлов. Для удобства работы названия признаков заменены на короткие. Итогововый датасет можно скачать <a href =\"https://drive.google.com/open?id=1BMKpbFUymeW4f16vc-HoChNxHewjJcXm\">тут</a>.\n",
    "\n",
    "### Признаки ###\n",
    "#### Целевой ####\n",
    "- __visitor:__ Число посещений – всего, измеряется в тыс.чел  \n",
    "\n",
    "#### Объясняющие ####\n",
    "- __Name:__ Название музея  \n",
    "- __kopuk:__ КОПУК (код классификатора по учреждениям культуры)\n",
    "- __adres:__ Адрес  \n",
    "- __pred_general:__ Число предметов основного фонда на конец года, единиц  \n",
    "- __pred_second:__ Число предметов научно-вспомогательного фонда на конец года, единиц  \n",
    "- __pred_elcat:__ Число музейных предметов, внесенных в электронный каталог музея, единиц  \n",
    "- __ex:__ Число выставок – всего, единиц  \n",
    "- __workers:__ Численность работников - всего, человек  \n",
    "- __build:__ Число строений - всего, единиц  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Соберем DataFrame из нескольких csv файлов. \n",
    "for i in range(9):\n",
    "    data = pd.read_csv(\"../../data/museum_research/data-\" + str(i+1) + \"-structure-1.csv\",)\n",
    "    data['year'] = 2008 + i #добавим еще признак, отражающий ГОД, на конец которого представлены данные\n",
    "    if i == 0:\n",
    "        df = data\n",
    "    else:  \n",
    "        df = pd.concat([df,data],axis=0, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('DataFrame содержит ' + str(df.shape[0]) +' объектов и ' + str(df.shape[1]) + ' признаков. Взглянем на них:')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В глаза бросается адрес Музея, как кладязь информации. \n",
    "Попробуем сразу получить больше информации из признака \"Адрес\", сгенерируем дополнительно признаки из географических координат:  \n",
    "- __longitude__: долгота\n",
    "- __latitude__: широта  \n",
    "Попутно из геоданных заберем пару категориальных признаков:\n",
    "- __region__: Название Субъекта\n",
    "- __okrug__: Название округа\n",
    "  \n",
    "Для геокодирования призанов, воспользуемся классом geocode из библиотеки __ArcGis__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gis = GIS()\n",
    "\n",
    "uni_adress = pd.unique(df['adres']) #выберем уникальные Адреса\n",
    "lat = {} #сформируем словарь для широты\n",
    "long = {} #сформируем словарь для долготы\n",
    "region = {}\n",
    "okrug = {}\n",
    "\n",
    "#Пробежимся по нашему списку уникальных адресов\n",
    "for i in tqdm_notebook(range(uni_adress.shape[0])): \n",
    "    geo = arcgis.geocoding.geocode(uni_adress[i], max_locations=1)[0] #получаем ГЕО данные \n",
    "    lat[str(uni_adress[i])] = geo['location']['y'] #формируем словари с координатами (x - долгота, y - широта)\n",
    "    long[str(uni_adress[i])] = geo['location']['x']\n",
    "    region[str(uni_adress[i])] = geo['attributes']['Region']\n",
    "    okrug[str(uni_adress[i])] = geo['attributes']['Territory']\n",
    "    \n",
    "df['latitude'] = df['adres']\n",
    "df['longitude'] = df['adres']\n",
    "df['region'] = df['adres']\n",
    "df['okrug'] = df['adres']\n",
    "\n",
    "df['latitude'] = df['latitude'].apply(lambda x: lat[x])\n",
    "df['longitude'] = df['longitude'].apply(lambda x: long[x])\n",
    "df['region'] = df['region'].apply(lambda x: region[x])\n",
    "df['okrug'] = df['okrug'].apply(lambda x: okrug[x])\n",
    "\n",
    "df.to_csv('../../data/museum_research/data-all-coord.csv',index_label='index') #сохраним данные"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> _Вычисления выше, занимают от 30 минут до часа, чтобы не собирать геоданные заного, данные сохранены в файл 'data-all-coord.csv'_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Пришло время, разделить данные на обучающую и контрольную выборки ##\n",
    "\n",
    "Напомню, что прогнозировать планируем посещения в 2016-ом году.\n",
    ">_Но сначала загрузим данные заного, если выше, что-то пошло не так :)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('../../data/museum_research/data-all-coord.csv',\n",
    "                 index_col='index') \n",
    "\n",
    "df_pred = df[df.year == 2016]\n",
    "df = df[df.year != 2016]\n",
    "print('Объектов для исследования: ', df.shape[0],', для прогноза:', df_pred.shape[0])\n",
    "print('Итоговый датасет содержит: ',df.shape[1], ' признаков, в том числе 1 целевой')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проводить исследование данных будем, принимая тот факт, что мы вообще не знаем данных за 2016 год, а лишь будем смотреть на нее как на отложенную выборку. \n",
    "\n",
    "##  Часть 2. Первичный анализ данных ##\n",
    "Разобьем признаки на типы и поговорим про каждый из них\n",
    "- Целевой \n",
    "    - visitor\n",
    "- Категориальные признаки: \n",
    "    - kopuk\n",
    "    - region \n",
    "    - okrug\n",
    "    - year\n",
    "- Вещественные признаки: \n",
    "    - pred_general \n",
    "    - pred_second\n",
    "    - pred_elcat\n",
    "    - ex\n",
    "    - workers\n",
    "    - build\n",
    "    - latitude\n",
    "    - longitude\n",
    "- Текстовые признаки: \n",
    "    - Name\n",
    "    - adres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Категориальные признаки: \n",
    "cat = ['kopuk','region','okrug','year']\n",
    "#Вещественные признаки:\n",
    "num = ['pred_general', 'pred_second', 'pred_elcat', 'ex', 'workers',\n",
    "       'build', 'latitude', 'longitude']\n",
    "#Текстовые признаки:\n",
    "txt = ['Name','adres']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Изучение целевого признака ###\n",
    "Целевая переменая, вещественное число и наша задача будет сведена к задаче восстановления регрессии, поэтому нам важно, чтобы целевая переменная имела нормальное распределение. \n",
    "\n",
    "Построим графики плотности рапределения целевой пременной, а так же сравним распредление с нормальным. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(rc={\"figure.figsize\": (20, 5)})\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "p_x = sns.distplot(df['visitor'], fit=st.norm, kde=True,ax=ax1, bins=50)\n",
    "ax1.legend()\n",
    "ax1.set_title('Плотность распределения целевой переменной')\n",
    "ax2 = fig.add_subplot(122)\n",
    "prob = st.probplot(df['visitor'], dist=st.norm, plot=ax2)\n",
    "ax2.set_title('Probabilyty plot')\n",
    "ax2.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из графиков становится видно, что распределение не симметрично, скошено и совсем не походит на нормальное распределение, Дополнительно проверим это статистическими тестами с помощью класса stats из библиотеки scipy, но для начала вспомним теорию.\n",
    "> _Распределение оценивается как предположительно близкое к нормальному, если установлено, что от 50 до 80 % всех значений располагаются в пределах одного стандартного отклонения от среднего арифметического, и коэффициент эксцесса по абсолютной величине не превышает значения равного двум._\n",
    "\n",
    "> _Распределение считается достоверно нормальным если абсолютная величина показателей асимметрии и эксцесса меньше их ошибок репрезентативности в 3 и более раз._\n",
    "\n",
    "> _Дополнительно применим тест Шапиро-Уилка_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_testing(vector):\n",
    "    testA = np.sqrt(6/len(vector))\n",
    "    testE = 2 * np.sqrt(6/len(vector))\n",
    "    A = np.abs(np.round(st.skew(vector),decimals=2))\n",
    "    E = np.abs(np.round(st.kurtosis(vector) ,decimals=2))\n",
    "    print('Асимметрия: ', A, '. Скошенность больше ошибки репрезентативности в : ',\n",
    "        np.abs(np.round(st.kurtosis(vector)/testA, decimals=2)) , ' раз.')\n",
    "    print('Эксцесс: ', E, '. Куртозис больше ошибки репрезентативности в : ',\n",
    "          np.abs(np.round(st.kurtosis(vector)/testE, decimals=2)) , ' раз.')\n",
    "    print('Гипотеза о нормальности распределения не отвергается с вероятностью: ', st.shapiro(vector)[1]*100, '%')\n",
    "    return \n",
    "\n",
    "norm_testing(df['visitor'].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Вывод 1: ###\n",
    "__Распределение не нормально__, повторюсь, что перед нами стоит задача восстановения регресии, то нам очень важно приблизить распределение целевой переменной к нормальному. (См. часть \"Предобработка данных\")\n",
    "\n",
    "### Изучение категориальных признаков ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[cat].head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(cat)):\n",
    "    print('Признак', cat[i], ' содержит ', len(np.unique(df[cat[i]].fillna(\"\"))), ' уникальных значений')\n",
    "    \n",
    "#fillna применён, потому что в признаке okrug есть какое-то пустое значение,\n",
    "# вызваное следствием ошибки в геокоддировании (См. часть \"Предобработка данных\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- признак kopuk, исходя из своей природы, является уникальным ID для каждого музея и имеет место  OneHotEncoding\n",
    "- признаки region и okrug заданы текстово и перед применением OneHotEncoding необходимо произвести label-кодирование\n",
    "- Признак year в рамках проекта будем рассматривать как вещественный, тем более, что в на отложенной выборке знаечение 2016 появляется впервые."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Категориальные признаки: \n",
    "cat = ['kopuk','region','okrug']\n",
    "#Вещественные признаки:\n",
    "num = ['pred_general', 'pred_second', 'pred_elcat', 'ex', 'workers',\n",
    "       'build', 'latitude', 'longitude','year']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Вывод 2: ###\n",
    "__Проведем OneHotEncoding по всем категориальным признакам, кроме year __\n",
    "\n",
    "### Изучение вещественных признаков ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[num].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Вывод 3: ###\n",
    "Вещественные признаки имеют разную размерность, поэтому необходимо произвести их стандартизацию \n",
    "### Изучение Текстовых признаков ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[txt].head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Название музея является таким же уникальным ID для каждого музея как kopuk, и будет с ним корреллировать,\n",
    "- Из признака adres, мы получили достаточно информации на этапе сбора датасета,\n",
    "В рамках текущего проекта предлагается данне признак не использовать.\n",
    "\n",
    "Однако, в дальнейшем их можно использовать в качестве источника для порождения дополнительных признаков (см. часть \"Выводы\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 3. Первичный визуальный анализ данных ##\n",
    "\n",
    "В рамках проекта решается задача прогнозирования посещений в 2016-ом году, посмотрим, на тенденцию изменения посещений по Годам."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pvt = df.pivot_table(columns=['year'], values='visitor', aggfunc='sum')\n",
    "plt.figure()\n",
    "sns.barplot(pvt.columns.values, pvt.values[0], color='Orange', palette=None,)\n",
    "plt.xlabel('Год')\n",
    "plt.ylabel('Кол-во посещений, тыс. чел.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Взглянем, на матрицу корреляции признаков:')\n",
    "sns.set(palette='GnBu_r', rc={\"figure.figsize\": (20, 5)})\n",
    "sns.heatmap(df.corr(), vmin=0, vmax=1,annot=True, cmap='Greens',)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Выводы: ###\n",
    "- Признак pred_elcat имеет сильную линейную связь (r=0.82) с pred_general, что вполне логично, исходя из природы данных. Возникает риск мультиколлениарности. \n",
    "- координаты музея, слабо коррелируют с целевым признаком. \n",
    "- целевой признка имеет наибольшую линейную связь (r=0.73) с workers, однако здесь нужно учесть, что скорее всего именно кол-во посещений влияло на кол-во сотрудников. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Стоит сделать предположение, что координаты, в некоторых строках определены не верно и данные являются выбросами\n",
    "\n",
    "Попробуем построить HeatMap посещений на карте РФ по координатам с помощью модуля Folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tunguzka = [60.7308, 101.9744]\n",
    "folium.Map.choropleth\n",
    "h_map = folium.Map(location=tunguzka, zoom_start=3,tiles='OpenStreetMap', width='100%', height = 700)\n",
    "\n",
    "data = [[x[0], x[1],x[2]] for x in np.array(df[df['year'] ==2015][['latitude','longitude','visitor']])]\n",
    "\n",
    "plugins.HeatMap(data, radius = 12,blur = 9).add_to(h_map)\n",
    "h_map.render()\n",
    "#h_map.save('map.html')\n",
    "h_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Судя по карте, можно предположить, что кол-во посещений сильно зависит от расстояния до центральных городов, как например Москва и Петербург и может послужить дополнительным признаком (см. часть \"Создание новых признаков и описание этого процесса\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 4. Закономерности, \"инсайты\", особенности данных##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Выбросы, пропуски, ошибки Гео-данных###\n",
    "Посмотрим, что из себя представляют выбросы на координатах:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Всего объектов:', df[df.latitude <= 25]['region'].shape[0],\n",
    "      ' в регионах:', np.unique(np.unique(df[df.latitude <= 25]['region'])))\n",
    "\n",
    "print('Всего объектов:', df[df.longitude <= 0]['region'].shape[0],\n",
    "      ' в регионах:', np.unique(np.unique(df[df.longitude <= 0]['region'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для Чукотского автономного округа такие координаты являются вполне нормальными, однако, Сингапура в РФ точно нет, судя по Wiki :)  \n",
    "\n",
    "Посмотрим, что это за музей втесался в наши официальные данные МинКульта:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.unique(df[df.region == 'Singapore']['Name'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ошибка в гео данных \"выскочила\" по одному музею, и можно без особого труда восстановить данные (см. часть \"Предобработка данных\")\n",
    "Правильные значения (ист. Wiki):\n",
    "- okrug: 'Центральный федеральный округ'\n",
    "- latitude: 55.7744\n",
    "- longitude: 37.6395\n",
    "- region: 'Москва'   \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### Проверка на выбросы по целевому признаку ###\n",
    "Проверим, что за аномальные скачки были у целевого признака и являются ли они выбросами"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "un_name = np.unique(df[df.visitor > 3000]['Name']) \n",
    "for i in range(len(un_name)):\n",
    "    print(un_name[i], ' за год посещают в среднем: ', \n",
    "          np.round(np.mean(df[df.Name == un_name[i]]['visitor']), decimals=2)\n",
    "          , 'тыс., а максимум достигал: ',\n",
    "            np.round(np.max(df[df.Name == un_name[i]]['visitor']), decimals=2), 'тыс.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Так как в выборку попали самые крупные музеи РФ, то смело можно сделать вывод, что данные скачки обусловленны природой данных, а не аномальными выбросами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 5. Выбор метрики##\n",
    "\n",
    "В задачах регрессии обычно используется среднеквадратичная ошибка или средняя абсолютная ошибка, которые подходят для сравнения двух моделей или для контроля качества во время обучения, но не позволяет сделать выводов о том, насколько хорошо данная модель решает задачу.\n",
    "Учитывая широкий разброс целевой переменной, и специфику задачи логично в качестве метрики выбрать коэффициент детерминации (r2_score) и посмотреть насколько хорошо модель объясняет данные или же наш прогноз сопоставим по качеству с константным предсказанием."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 6. Выбор модели ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Предлагается использовать линейную модель для восстановления регрессии.  \n",
    "\n",
    "__Предпосылки: __\n",
    "- Большое кол-во признаков, с разряженными матрицам\n",
    "- Хочется интерпретировать коэфициенты для выявления влияния факторов\n",
    "\n",
    "__Ограничения: __\n",
    "- Имеем риск получения плохого результата, так как зависимость целевой переменной от признаков нелинейная (Видно из матрицы корреляции)\n",
    "\n",
    "В качестве линейной модели будет использован __SGDRegressor__:\n",
    "- Быстро обучается\n",
    "- Гибко настраивается\n",
    "- Позволяет использовать регуляризацию"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 7. Предобработка данных ##\n",
    "\n",
    "\n",
    "### Исправление ошибок Геокодинга ###\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[df.region == 'Singapore']['okrug'] = 'Центральный федеральный округ'\n",
    "df[df.region == 'Singapore']['latitude'] = 55.7744\n",
    "df[df.region == 'Singapore']['longitude'] = 37.6395\n",
    "df[df.region == 'Singapore']['region'] = 'Москва'\n",
    "\n",
    "df_pred[df_pred.region == 'Singapore']['okrug'] = 'Центральный федеральный округ'\n",
    "df_pred[df_pred.region == 'Singapore']['latitude'] = 55.7744\n",
    "df_pred[df_pred.region == 'Singapore']['longitude'] = 37.6395\n",
    "df_pred[df_pred.region == 'Singapore']['region'] = 'Москва'\n",
    "\n",
    "# !!!!!!!!!!!! Не присваиваются значения! Исправь! \n",
    "df.dropna(inplace=True)\n",
    "df_pred.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Обработка целевого признака для приведения к нормальному распределению ###\n",
    "Попробуем для начала просто взять логарифм из целевой переменной"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_modified = np.log((df['visitor'].values)+1)\n",
    "norm_testing(target_modified) #Используем тест, который использовали ранее\n",
    "\n",
    "#Построим графики\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "prob = sns.distplot(target_modified, fit=st.norm, kde=True,ax=ax1)\n",
    "ax1.set_title('Гистограмма распределения целевой переменной')\n",
    "ax2 = fig.add_subplot(122)\n",
    "prob = st.probplot(target_modified, dist=st.norm, plot=ax2)\n",
    "ax2.set_title('Probabilyty plot')\n",
    "ax2.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Логарифмирование дало эффект, однако, для достижения большего результата применим Метод Бокса-Кокса, так же из класса stats из библиотеки scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "step_boxcox = 0.592\n",
    "target_modified = (df['visitor']+step_boxcox)\n",
    "par = st.boxcox(target_modified)\n",
    "target_modified = par[0]\n",
    "par_lamb = par[1] #Пригодится, для изменения тестовой выборки\n",
    "norm_testing(target_modified) #Используем тест, который использовали ранее\n",
    "df['target'] = target_modified\n",
    "#Построим графики\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "prob = sns.distplot(target_modified, fit=st.norm, kde=True,ax=ax1)\n",
    "ax1.set_title('Гистограмма распределения целевой переменной')\n",
    "ax2 = fig.add_subplot(122)\n",
    "prob = st.probplot(target_modified, dist=st.norm, plot=ax2)\n",
    "ax2.set_title('Probabilyty plot')\n",
    "ax2.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Картина стала в целом лучше, но Говорить о нормальности распредления нам не позволяет тест Шапиро-Уилка.\n",
    "> _Прим.: в документации scipy к данному тесту, есть ссылка, что при более чем 5000 объектов p-value может отображаться не корректно._\n",
    "\n",
    "Применим Трансформацию и к отложенной выборке:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pred['target'] = df_pred['visitor']+step_boxcox\n",
    "df_pred['target'] = st.boxcox(df_pred['target'].values, lmbda=par_lamb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Удалим старую целевую переменную\n",
    "df.drop('visitor', inplace=True, axis=1)\n",
    "df_pred.drop('visitor', inplace=True, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Проведем OneHotEncoding по всем категориальным признакам ###\n",
    "\n",
    "Так как уникальные значения категориатных признаков в 2016-ом году встречаются не все, то чтобы ничего не потерять, соединим наши выборки, а потом снова разделим."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_u = pd.concat([df, df_pred])\n",
    "df_OHE = pd.concat([df_u.drop('kopuk', 1), pd.get_dummies(df_u['kopuk'],prefix='kopuk') ], axis=1).reindex()\n",
    "df_OHE = pd.concat([df_OHE.drop('region', 1), pd.get_dummies(df_OHE['region'],prefix='reg') ], axis=1).reindex()\n",
    "df_OHE = pd.concat([df_OHE.drop('okrug', 1), pd.get_dummies(df_OHE['okrug'],prefix='okr') ], axis=1).reindex()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Удалим Текстовые фичи ###        \n",
    "\n",
    "Обоснование удаления приведено в части \"Первичный анализ данных\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_txt = df_OHE.drop(['Name','adres'],axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Шкалирование вещественных признаков ###\n",
    "\n",
    "Шкалирование вещественных признаков и повторное разбиение на тестовую и отложеную выборки будет произведено в части \"Кросс-валидация, подбор параметров\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 8. Создание новых признаков и описание этого процесса ##\n",
    "\n",
    "Ранее выдвигалось предположение, о том, что расстояния до Москвы и Санкт-Петербурга могут влиять на кол-во посещений. За основу расчета растояний возьмем Координаты Московского кремля и Дворцовой площади"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kremlin = [55.750730, 37.617322] # Координаты Московского кремля\n",
    "dp_piter = [59.9384, 30.3227] #Координаты Дворцовой площади\n",
    "dist_to_kremlin = []\n",
    "dist_to_dp = []\n",
    "\n",
    "for index, row in tqdm_notebook(df_txt.iterrows()):\n",
    "    dist_to_kremlin = np.append(dist_to_kremlin, dist.great_circle(( row['latitude'], row['longitude']),(kremlin[0],kremlin[1])).km)\n",
    "    dist_to_dp = np.append(dist_to_dp, dist.great_circle(( row['latitude'], row['longitude']),(dp_piter[0],dp_piter[1])).km)\n",
    "    \n",
    "df_txt['dist_to_kremlin']=dist_to_kremlin\n",
    "df_txt['dist_to_dp']=dist_to_dp\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Рассмотрим облако рассеивания новых переменных и целевой"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.PairGrid(df_txt[['target', 'dist_to_kremlin','dist_to_dp']], size=5)\n",
    "g = g.map(plt.scatter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Из графика видно, что связь не линейна, чтобы в этом убедиться взглянем на матрицу корреляций:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_txt[['target', 'dist_to_kremlin','dist_to_dp']].corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Итак:__ Новые признаки не имеют линейной связи, но отрицательные коэфициенты корреляции говорят о наличии обратной зависимости: Чем больше расстояние, тем меньше целевая переменная, что соответсвует нашему изначальному предположению."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 0. Финальная предобработка данных ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Повторно разделим выборку на Обучение и Тест\n",
    "df_final_pred = df_txt[df_txt.year == 2016]\n",
    "df_final_train = df_txt[df_txt.year != 2016]\n",
    "print('Объектов для обучения: ', df_final_train.shape[0],', для прогноза:', df_final_pred.shape[0])\n",
    "print('Итоговый датасет содержит: ',df_final_train.shape[1], ' признаков, в том числе 1 целевой')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = df_final_train.drop('target', axis=1).values\n",
    "y_train = df_final_train['target'].values\n",
    "\n",
    "X_test = df_final_pred.drop('target', axis=1).values\n",
    "y_test = df_final_pred['target'].values\n",
    "\n",
    "#Проведем стандартизацию признаков, о котором говорили в части \"Предобработка данных\"\n",
    "scaling = sklearn.preprocessing.StandardScaler()\n",
    "X_train = scaling.fit_transform(X_train)\n",
    "X_test = scaling.transform(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 9. Кросс-валидация и настройка гиперпараметров модели ##\n",
    "\n",
    "Части \"Выбор модели\" было описано, что в качестве модели мы будем использовать стохастический градиентный спуск. \n",
    "Принимая во внимание, что данные имеют большие выбросы, я буду использовать функцию потерь Губера, её использование позволит уменьшить вклад выбросов и большой квадратичной ошибки, а неточности будут проигнорированы. \n",
    "\n",
    "\n",
    "Крос-валидация будет проведена на трех фолдах с перемешиванием.  \n",
    "Гиперпарметрами для функции будет: \n",
    "- Alpha - определяющий порог регуляризации\n",
    "- Epsilon - порог, при котором становится менее важным получить предсказание точно правильно. Тоесть когда функция потерь изменит свою работу с вида L<sub>2</sub> на L<sub>1</sub>\n",
    "\n",
    "Обучим класс GridSearchCV у четом сказанного выше."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kf = model_selection.KFold(n_splits=3, shuffle=True, random_state=18)\n",
    "\n",
    "SGD = linear_model.SGDRegressor(loss='huber', max_iter=100, penalty='l2',\n",
    "                                random_state=18,learning_rate='optimal')\n",
    "\n",
    "params = {'alpha': np.linspace(0.01,0.04, num=5),\n",
    "          'epsilon': np.linspace(0.1,2, num=5)}\n",
    "\n",
    "grid = model_selection.GridSearchCV(SGD, param_grid=params, scoring='r2', cv=kf, verbose=3)\n",
    "grid.fit(X_train,y_train)\n",
    "grid.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid.best_estimator_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На кросвалидации R<sup>2</sup> - составил 0.86, при параметрах:  \n",
    "- alpha = 0.04  \n",
    "- epsilon = 1.525\n",
    "\n",
    "Учитывая, что коэффициент детерминации выше 80%, можно говорить о хорошем качестве модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 10. Построение кривых валидации и обучения ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = np.linspace(0.01,0.04, num=5)\n",
    "val_train, val_test = model_selection.validation_curve(SGD, X_train, y_train, param_name='alpha',\n",
    "                                                       param_range=alphas, cv=kf,\n",
    "                                       scoring='r2')\n",
    "\n",
    "def plot_with_err(x, data, **kwargs):\n",
    "    mu, std = data.mean(1), data.std(1)\n",
    "    lines = plt.plot(x, mu, '-', **kwargs)\n",
    "    plt.fill_between(x, mu - std, mu + std, edgecolor='none',\n",
    "                     facecolor=lines[0].get_color(), alpha=0.2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_sizes = np.linspace(0.05, 1, 10)\n",
    "N_train, tr_train, tr_test = model_selection.learning_curve(SGD,\n",
    "                                                  X_train, y_train, train_sizes=train_sizes, cv=kf,\n",
    "                                                  scoring='r2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(20, 5))\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.set_title('Кривая Валидации')\n",
    "plot_with_err(alphas, val_train, label='training scores')\n",
    "plot_with_err(alphas, val_test, label='validation scores')\n",
    "plt.xlabel(r'$\\alpha$'); plt.ylabel('r^2')\n",
    "plt.legend();\n",
    "\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "ax2.set_title('Кривая Обучения')\n",
    "plot_with_err(N_train, tr_train, label='training scores')\n",
    "plot_with_err(N_train, tr_test, label='validation scores')\n",
    "plt.xlabel('Training Set Size'); plt.ylabel('r^2')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Валидационная кривая:\n",
    "- Так кривые валидации близко друг к другу и их значение близко к 1, то можно говорить об отсутствии переобучения или недообучения. \n",
    "\n",
    "Кривая обучения:\n",
    "- Так как кривые сошлись друг к другу, добавление новых данных ничего не изменит."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Часть 11. Прогноз для тестовой выборки ##\n",
    "\n",
    "Спрогнозируем колличество посещений в музея в 2016 году по результатам модели и расчитаем кофициент детерминации для Тестовой выборки:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(metrics.r2_score(y_test, grid.predict(X_test)), decimals=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "R<sub>2</sub> Для тестовой выборки соотносится с показателем достигнутым на кросс-валидации. \n",
    "Можно сделать вывод о том, что модель имеет хорошее качество.  \n",
    "Посмотрим, какие факторы максимально влияют на целевую переменую:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_coef = pd.DataFrame({ 'name': df_final_train.drop('target', axis=1).columns, 'value': grid.best_estimator_.coef_ })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_top_features = 15\n",
    "coef_plot = pd.concat([df_coef.sort_values(by='value', ascending=False).head(n_top_features),\n",
    "                        df_coef.sort_values(by='value', ascending=False).tail(n_top_features)], axis=0)\n",
    "\n",
    "interesting_coefficients = coef_plot['value']\n",
    "feature_names = coef_plot['name']\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "colors = [\"red\" if c < 0 else \"blue\" for c in interesting_coefficients]\n",
    "plt.bar(np.arange(1, 2 * n_top_features+1), interesting_coefficients, color=colors)\n",
    "plt.xticks(np.arange(1, 2 * n_top_features+1), feature_names, rotation=90, ha=\"center\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Часть 12. Выводы ##\n",
    "\n",
    "Создана модель, которая имеет хорошую обобщающую способность. \n",
    "Ценность модели заключается в том, что ее можно использовать для прогнозирования кол-ва посетителей в музеях. На основании полученных данных можно оптимизировать затраты. \n",
    "\n",
    "Кроме того, можно выделять факторы, которые максимально влияют на посещаемость музея\n",
    "\n",
    "В части путей улучшения, можно использовать следующие направления:\n",
    "- Порождение дополнительных фич из названия музея, например его категорию (Археологический, Исторический, Краеведческий) или оценку (отзывы поситителей).\n",
    "- Порождатьть дополнительные фичи на основании геоданных из Адреса.\n",
    "- Для повышения интерпритируемости, возможно Смещение фич на +1 Год для прогноза, так как датасет содержить данные на конец года, но формально они соответствуют данным на начало года"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
